
function [ES1,EQ1,S1,EQU1] = surplus1(gamma,B,D,P,xi,n,options,cmax,save_eqL,save_eqU)

N = sum(n);
NS = size(xi,1);

% PRI-PVEM joint surplus with joint PRI candidate
ES1 = zeros(N,NS);
if save_eqL == 1
    EQ1 = zeros(N,5,NS);
    S1 = zeros(N,5,NS);
end
if save_eqU == 1
    EQU1 = zeros(N,5,NS);
end

for t = 1:N
    Xt = [blkdiag(kron(eye(2),[0,1,0,D(t,:)]),kron(eye(2),[0,1,D(t,:)]),...
        [0,1,0,D(t,:)]),log([P(t,1:2),sum(P(t,3:4)),sum(P(t,3:4)),P(t,5)]')];
    parfor ns = 1:NS
        if save_eqU == 1
            [eqL,~,eqU,~] = EQ_bounds(gamma,B,Xt,xi(ns,:)',options,cmax,1); % largest and smallest spending equilibria
        else
            eqL = EQ_bounds(gamma,B,Xt,xi(ns,:)',options,cmax,1);
        end
        [~,s] = Shares([B;0;0;xi(ns,:)'],zeros(5,1),[Xt,eqL,eqL.^2],0,ones(5,5),n,zeros(1,2),1);
        ESt(ns) = (gamma(3) + gamma(4)) * log(s(4)) - eqL(4);
        if save_eqL == 1
            EQt(:,ns) = eqL;
            St(:,ns) = s;
        end
        if save_eqU == 1
            EQUt(:,ns) = eqU;
        end
    end
    ES1(t,:) = ESt;
    if save_eqL == 1
        EQ1(t,:,:) = EQt;
        S1(t,:,:) = St;
    end
    if save_eqU == 1
        EQU1(t,:,:) = EQUt;
    end
end



